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ABSTRACT 

We present a new and completely general technique for calculating the fine-grained phase- 
space structure of dark matter throughout the Galactic halo. Our goal is to understand this 
structure on the scales relevant for direct and indirect detection experiments. Our method is 
based on evaluating the geodesic deviation equation along the trajectories of individual DM 
particles. It requires no assumptions about the symmetry or stationarity of the halo formation 
process. In this paper we study general static potentials which exhibit more complex behaviour 
than the separable potentials studied previously. For ellipsoidal logarithmic potentials with 
a core, phase mixing is sensitive to the resonance structure, as indicated by the number of 
independent orbital frequencies. Regions of chaotic mixing can be identified by the very rapid 
decrease in the real space density of the associated dark matter streams. We also study the 
evolution of stream density in ellipsoidal NFW halos with radially varying isopotential shape, 
showing that if such a model is applied to the Galactic halo, at least 10^ streams are expected 
near the Sun. The most novel aspect of our approach is that general non-static systems can 
be studied through implementation in a cosmological N-body code. Such an implementation 
allows a robust and accurate evaluation of the enhancements in annihilation radiation due to 
fine-scale structure such as caustics. We embed the scheme in the current state-of-the-art code 
GADGET-3 and present tests which demonstrate that N-body discreteness effects can be kept 
under control in realistic configurations. 
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1 INTRODUCTION 

Dark Matter (DM) is still to be directly detected. The first indirect 
indications of its existence came in the 1930s, with measurements 
of the velocities of galaxies in clusters. The cluster mass required to 
gravitationally bind the galaxies was found to be more than an order 
of magnitude larg er than the sum of the lumin ous masses of the in- 
dividual galaxies jZwickvll 19331 : ISmithlll936h . The early detection 
of large amounts of uns een matter associated with the Local Group 
jKahn&Woitieill 19591) was followed in the 1970s by observations 
of the rotation curves of spiral galaxies which showed that these are 
flat, or eve n rising, at dist ances far beyond their stellar component s 
jRubin & Ford 1970; Fab er & Gallaghejl 19791 : iRubin et alj|l980h . 
Studies of satellite systems suggested that the mass distributions 
of most galaxies might be an or der of magnitude la r ger and more 
massi ve than their visible parts jOstriker et al.ll 19741 ; Einasto et al.l 
11974) . All these discoveries led to the conclusion that a large frac- 
tion of mass in the Universe is dark. This has also been supported 
by recent gravitational lensing studies that demonstrate the ex- 
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istence of extended massive DM halos (e.g. Mandelbaum et al. 
2006). 

As the dominant mass component of galaxies and large-scale 
structures, DM has necessarily become a key ingredient in theo- 
ries of cosmic structure formation. The most successful of these 
theories is the hierarchical paradigm. In the current version of this 
model, the DM is composed of nonbar yonic elementa ry particles 
known as Cold Dark Matter (CDM) Ipeeblesi [l983) . The term 
"cold" derives from the fact that the DM particles had negligi- 
ble thermal motions at the time of matter-radiation equality. Their 
abundance was set when the interaction rate became too small for 
them to remain in thermal equilibrium with other species in the 
expanding Universe. The first objects to form in a CDM Universe 
are small galaxies. They then merge and accrete to give rise to the 
larger structures we observe today. Thus structure formation oc- 
curs in a "bottom-up" fash ion tBlumenthal et al.lll984l ; lDavis et al.l 
ll985l ; ISpringeletalj2006l) . 

The crucial test of this paradigm undoubtedly consists in the 
determination of the nature of DM through direct detection exper- 
iments. Among the most promising candidates from the particle 
physics perspective are axions and neutralinos. Axions were intro- 
duced to explain the absence i n Nature of stroiig-CP (Charge con- 
jugation and Parity) violations l lPeccei 8i Ouin rJl977h . The cosmo- 
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logical population formed out of equilibrium as a zero momentum 
Bose condensate. They can be detected through their conversion to 
photons in the presence of a strong magnetic field, (e.g. Ogawa et 
al. 1996; Hagmann et al. 1998). Neutralinos are the lightest super- 
symmetric particles, and may be considered as a particular form of 
weak ly interacting massive particles (WIMPs) jSteigman & Turner! 
[l983). The most important direct detection process for neutralinos 
is through elastic scattering on nuclei. 

Today many exper iments ar e search i ng for these parti- 
cles. jAkerib&et all 12004 : Sanglar d & et alll2005l ; ISchnedliood : 
lAprile et al.l 20071 : ISDOonerii2007.) . The main challenge lies in the 
large background they encounter. Having an idea what the detector 
signal might look like can help substantially in fine-tuning the ex- 
periments in order to increase the chance of a detection. In addition, 
many experiments are attempting to detect WIMP s indirectly b; 



searching for 7-ray emission from their annihilation ( de Boer 200: 
de Boer et al. 2005; Bergstrom & Hooper 2006; Hooper & Serpico 
2007). Predictions for this radiation are currently uncertain because 
very substantial enhancements are possible, at least in principle, 
from fi ne-scale structure such as caustics in the dark ma tter distri- 
bution ('Palcanton & HogarJl200ll :lMohavaee e t al.ll2007l) . 

The differential detection rate for WIMPs in a given material 
(per unit detector mass) is sensitive to the local velocity distribu- 
tion: 
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where Q is the energy deposited in the detector and da / dQ is the 
differential cross section for WIMP elastic scattering with the target 
nucleus. is the WIMP mass which lies in the GeV to TeV range 
for currently preferred models, p is the local mass density and /(v) 
the velocity distribution of WIMPs that reach the detector. 

From Eq. ([TJ it is evident that the count rate in a direction- 
insensitive experiment depends on the velocity distribution of the 
incident particles and wil l be modulated by t he orbital motion of 
the Earth around the Sun i iDrukieretal.lll986l) . In most studies, an 
isotropic Maxwellian distribution relative to the galactic halo has 
been assumed, (e.g. Freese et al. 1988), although there are other 
mode ls in the literature, discussing, for example, multivariate Gaus- 
sians jEvans et alj200(lh . Some attempts at understanding the effect 
of fine-scale structur e in the velocity distribution of DM particle s 
have also been made ( ISikivielll998l : IStiff et alj200ll : lHoganll200ll) . 

A significant sig nal could come from what are known 
as s treams of DM (' Sikivie et al] 1 19951 : iHelmi et al.l |2002| : 
iNataraian & Sikivie 2005.) . Both axions and WIMPs are cold. In the 
absence of clustering their present day velocity dispersion would be 
neghgible (Sv ~ 10~^°c for WIMPs and (5u ~ 10~^^c for axions). 
They are effectively restricted to a 3D hypersurface, a sheet in 6D 
phase-space. The growth of structure results in continual stretching 
and folding in phase-space of this initially almost uniform sheet. 
This process is called mixing. The more strongly the system mixes, 
the more streams of DM particles are present at a given location 
in configuration-space. Mixing stretches each sheet and so its den- 
sity decreases. The existence of distinct streams is a direct con- 
sequence of the coUisionless character and the coldness of CDM. 
The situation is sketched in Fig.[T| At the points where the number 
of streams changes by two, the local configuration-space density 
of dark matter becomes extremely high. These are caustics of th e 
kind studied by catastrophe theory dOilmorel I98l : lTremainel 199^ . 
Note that Liouville's Theorem prevents the CDM phase-space sheet 
from ever tearing, although it can be arbitrarily strongly stretched. 

If the dark matter density in the solar neighbourhood is dom- 




Figure 1. Sketch of an idealised CDM phase-space sheet in the x, x plane. 
The thickness of the line represents the local velocity dispersion within each 
stream. The small wiggles correspond to initial density perturbations and 
the multi-valued region reflects the multiple streams created by winding in 
non-linear regions. Depending on x position an observer sees one or three 
streams. At points where the number of streams changes by two, a caustic 
with a very high DM density is present. Such caustics may be significant for 
the total annihilation fiux. The number of streams at each point is a measure 
of the local amount of mixing. The dashed line represents the Hubble Flow. 
The cross marks the phase-space coordinates of a particular CDM particle 
embedded in the flow. 



inated by a small number of streams, its velocity distribution will 
effectively consist of a few discrete values, one for each stream. 
If, on the other hand, the number of streams is very large, the ve- 
locity distribution will be smooth and individual streams will be 
undetectable. This issue has so far been addressed only under sim- 
plifie d conditions and divorced from its proper c osmological con- 
text ( Helmi et al.ll2002l : INataraian & Sikiviell2005il . This is largely 
because the only tool capable of studying cosmological structure 
formation in full generality, namely N-body simulations, cannot re- 
solve the relevant scales. For example, Natarajan & Sikivie ( 20qJ) 
estimate that of the order of lO^'^ particles would be required in a 
simulation of the Milky Way's halo to begin to resolve the streams 
in the solar neighbourhood. E ven the largest simula tion so far pub- 
lished, the Via Lactea model ( Diema nd et 'Zl l2007h . is 4 orders of 
magnitude short of this minimum requirement. It is thus impossible 
to figure out the number of streams near the Sun and the properties 
of "typical" streams with current N-body capabilities. 

A related issue that has been much discussed is whether 
caustics can affect direct or indirect detection experiments. Some 
authors claim that caustics play an important role as stable 
phenomena connected to any collapsing CDM h alo and taking 
the form o f relatively massive rings or shells jSikivid [T999t 
lNataraian & Sikivie 2006). Other authors have argued that caus- 
tics should be mor e abundant, weaker and dynamically negligible 
jHelmi et alj|2002l) . Even if caustics negligibly affect the gravi- 
tational potential, they may have very substant i al effects on the 
annihil ation rate of dark matte r jHogan 2001; Bergstrom et ajj 
200 1I: IPieri & Branchinil l2005l: iMo havaee & Shandarin 200i 
Mohavaee et al.ll2007l : Nataraiar 2007b . All these papers were able 



to evaluate the enhancements due to caustics only under restrictive 
and unrealistic assumptions about symmetry, formation history, etc. 
While they demonstrate that large enhancement factors may be pos- 
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sible, they do not provide reliable estimates for the actual enhance- 
ment expected. The method we present below is capable of provid- 
ing such estimates for halos growing as predicted by the ACDM 
model. 

In this paper we will present a novel approach that directly 
analyses structure in the fine-grained phase-space distribution. 
We circumvent the "number of particle" problem by solving the 
geodesic deviation equation (GDE) for every DM particle. This al- 
lows us to calculate the local properties of the DM stream each 
simulation particle is embedded in, in particular, its configuration- 
space density and its local velocity distribution. The mass-weighted 
number of streams near any point is then the total local density di- 
vided by the mean density of the individual local streams. Caus- 
tic passages can be detected robustly from the properties of each 
particle's local stream, and the particle's instantaneous annihilation 
probability within this stream can be evaluated simply and inte- 
grated accurately through caustics. 

This paper is the first of a series. Here we introduce our 
method and apply it to study the evolution of streams in static po- 
tentials that are too complex to be analysed using previous tech- 
niques. We also implement our scheme in a state-of-the-art N-body 
simulation code, and use simple test problems to demonstrate that 
N-body discreteness effects can be kept under control in realis- 
tic applications. Later papers will address issues associated with 
mixing, caustics and annihilation radiation in the full cosmological 
context. 

The outline of our paper is the following: In Section 2 we 
present a detailed derivation of the GDE and show how it can 
be used to quantify mixing, to locate caustics, and to calculate 
stream densities and annihilation rates. Section 3 describes our 
code, DaMaFlow, which is designed to solve the GDE for single 
orbits in a wide variety of potentials. In Section 4 we analyse static, 
separable potentials and compare results from our method to pre- 
vious work. The fifth section applies our scheme to non-integrable, 
but still static potentials, revealing their complex phase-space struc- 
ture. In Section 6 we turn to more realistic non-spherical CDM 
halo potentials and discuss the influence of halo shape on stream 
density behaviour. Section 7 discusses the implementation of our 
method in an N-body code and presents results of simple tests of 
when discreteness effects compromise studies of stream densities 
and caustics. The final section summarises our results and gives 
some conclusions. 



2 THE GEODESIC DEVIATION EQUATION 

Our scheme for calculating the evolution of the fine-grained phase- 
space distribution in the neighbourhood of a DM particle is based 
on the evolution of the distance between two infinitesimally sep- 
arated particle trajectories. This evolution is described by the 
Geodesic Deviation Equation (GDE). We use the following no- 
tation to clearly distinguish between three- and six-dimensional 
quantities: an underline denotes a 'M? vector and two of them de- 
note a R'^^'^ matrix. An overline denotes a R" vector and two of 
them denote a R''^'' matrix. Thus a general phase-space vector is 
composed of two three-dimensional vectors: x — {x, v_). 

To derive the GDE we first write down the equations of motion 
for a DM particle. These are simply 

£(i;£o.^o) = -^x*(£(i;£o.l^o)). (2) 
with initial conditions a; (to ; SIq , ) = Xq w (t o ; ' — o ) ~ — o • 



As x = v_ the equation of motion in phase-space can be writ- 



ten: 



rr(t;xo)= j^^{x{t;xo)), (3) 

with initial conditions x (to; ^o) —xo = (So' — o)- 

We want to calculate the local stream density around the DM 
particle whose trajectory in phase-space is given by x {t; xq). To do 
so, we first ask how the displacement vector to a neighbouring DM 
particle in phase-space evolves with time: 



5(t) = X [t; Xo + So) - X {t; Xq) . 



(4) 



Note that 5 (t) is the displacement between the reference DM par- 
ticle, which was at xq at time to, and another particle which was 
at Xo + So at to- We are interested in properties in the immediate 
neighbourhood of the reference particle, so So is an infinitesimal 
displacement. This allows us to work to linear order: 



5(t)e^ (5o- V:,o)x(t;xo). 



(5) 



Introducing the phase-space distortion tensor D (note that this is a 
6x6 tensor) 



D {t;xo) = {t;xo) 

OXo 



(6) 



we can rewrite Eq. l[5ll as a simple linear transformation from the 
starting phase-space displacement So to the displacement S{t) at 
time t: 



S {t) ^ D {t-xo) So. 



(V) 



Because So is an arbitrary displacement vector, the distortion ten- 
sor describes how the complete local phase-space neighbourhood 
around the reference DM particle gets distorted while it is orbiting 
in the given potential. The time evolution of S{t) follows from the 
time evolution of the two trajectories. Again we can work this out 
in linear orde^^: 

'B{t;xo)So = t(t) 

= '^{x{t;xo)+S{t)) -'^{x{t;xo)) 
= (S{t)-V,)^{x{t-xo)) 

[^{t;xo)5o) ■V^^{x{t;xo)). (8) 



To derive the equation of motion for the distortion tensor itself we 
evaluate Eq. ^ for six unit vector phase-space displacements Sq'' 
with S(/l^ = Saj where a, j = 1,2, ... ,6, and Sab is the Kronecker 
delta. Taking into account Einstein's sum convention this yields for 
each component of Eq. l[8j 

D,j{t;xo) = D^o. {t;xo) S^'l 

D^-, (t;xo)S^,'^l^^ ^,{x{t;xo)) 

= [Dfj^{t;xo)5^j^-J-^i{x{t;xo)) 



D. 



Hi 



{t;xo) 7^ ) {x{t;xo)) 



dxfi J 



= Tif3{t;xo) DfSj {t;xo) , 



^ As (5o can be chosen arbitrarily small, this is always possible. 



(9) 
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where we have introduced the phase-space tidal tensor 



and 



T{t-xo) = 







1 



(10) 



and T is the configuration-space tidal tensor given by the second 
derivatives of the gravitational potential Tij — —d^^/dxidxj. As 
we are are only interested in linear order we replace = by = in Eq. 
([9} and get an equation of motion for the 6D distortion tensor: 



D{t-xo) = T (t-xo) D {t;xo) . 



(11) 



Note that this first-order tensor differential equation represents a 
system of 36 coupled ordinary first-order differential equations. To 
solve them we need to specify initial conditions. These follow from 
the constraint 



So ^ D (to; a;o) So 



D (to;xo) 



(12) 



DM behaves like a coUisionless fluid and its fine-grained phase- 
space density f{x, t) is described by the well-known Vlasov equa- 
tion: 



9/ 
dt 



(13) 



Using the Lagrangian derivative this reads D//Dt = 0, meaning 
that the local fine-grained phase-space density has to be conserved 
along the orbit of every particle in the system. Imagine A*' particles 
that fill a small phase-space volume dVb = clxgdi^g around the 
DM reference particle at time to- At a later time t these particles 
fill a volume AV = AxAv^ around the reference particle at x{t\ xo)- 
Conservation of phase-space density implies that the two volumes 
are the same dVo = AV. The evolution of the initial displacement 

— (n) 

vectors Sq from the reference particle to one of the N other par- 
ticles in the volume is described by the distortion tensor associated 
with the reference particle: 



5'"' {t)^D{t-xo)S^o ^ 



(14) 



The change in volume due to this linear transformation is given 
by the determinant det(D {t;xo)). As this volume has to be con- 
served, the determinant of the phase-space distortion tensor has to 
be conserved. Note that not only must the volume be conserved, 
but also the orientation of the volume element. This means that the 
sign of the determinant is also fixed. From the initial conditions Eq. 
l ll2t one gets det(_D {t;xo))=i at all times. This fact can be used 
to check numerical solutions of the equations. 

The structure of Eq. dllb allows the equations of motion for 
the distortion to be broken down to a set of equations that is more 
convenient to work with. Let us first rewrite Eq. Jilt using blocks 
of 3 X 3 tensors Q 



D D 
At\D D 



D D 

^=XV ^=X 

D D 



g i 

T 2 

D 

VX _ _ 

T D T D 



D 



Writing down the equation for each matrix block yields four equa- 
tions: 



D 



D 



D 



D 



(15) 



D =T D 
These can be combined to give: 
b D 

X = X 



D =T D 



D =r D . 



(16) 



(17) 



Thus we get two identical differential equations of second-order 
for two 3x3 tensors whose dynamics is driven by the ordinary 
tidal tensor. From the initial condition for the 6D distortion one 
can see that the only difference between D and D lies in the 

appropriate initial conditions: -D (to;aio)~i'— (*o;2;q) = 
and_D {to', Ho) ~ 0, D_ (to; £o) = 1 in the two cases. From the 
solutions of these two initial condition problems the 6D distortion 
solution can then be constructed: 



D - 



D 



D 



A/ AW A/AtD 



(18) 



Up to this point we have worked out all equations in phase- 
space. As we are interested in the stream density in configuration- 
space, we need to project down to this space. As already mentioned, 
CDM lies on a thin sheet in phase-space. This sheet has a certain 
orientation at the starting point of the reference DM particle. Take 
{x,v) : V — V_{x;to,xo) to be the local parametrisation of the 
sheet surrounding this particle at time tfl Now we ask how an 
infinitesimal displacement So^x configuration-space is distorted 
by evolution. First we note that any displacement in x implies a 
displacement in velocity-space due to the restriction of particles to 
the sheet: 

_ _ dV _ 

io,« (a;o)io,^ ; ^^i^o) = ^ {xo-to,xo) . (19) 

The phase-space distortion D describes how the corresponding 
phase-space displacement [Sg ^,5g ^) evolves. We are here inter- 
ested in the configuration-space part of the phase-space displace- 
ment at time t: 

W (*;^o)^o.. (i;^o)Z, (^o)io,.- (20) 

The evolution of the displacement in configuration-space can also 
be described by a linear transformation: 

L (t)=mt;xo)S,-,,^, 

where we have introduced the configuration-space distortion tensor 
(note that this is a 3 x 3 tensor) 

D{t-xo) ^D^Jt-xo)+D^Jt-xo)^^{xo). (21) 

This tensor can also be derived with the help of two projection op- 
erators: 



We suppress the argument t; xq to avoid confusion. 



^{t-xo) = [1 g)D {t-xo) y |_ ^ j . (22) 

As in the case of phase-space distortion, the change in volume due 
to the linear transformation in Eq. | |2U is given by the determinant. 



Such a parametrisation is always possible locally, but due to mixing there 
is, in general, no simple global parametrisation of the stream. This is only 
possible for very early times, where the sheet is dominated by the Hubble 
Flow, the X — V relation is one-to-one, and the stream density is almost 
uniform. 
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Distortion tensor 
8(t) = P(t)8, 




Figure 2. The configuration-space distortion tensor 13 describes liow an 
initial small configuration-space displacement 5q evolves in time. This re- 
flects the stretching of an infinitesimally small cloud of virtual particles 
around the reference particle that is placed at Xq at time Iq . The stretching 
of the cloud is driven by the tidal field that the DM particle encounters as it 
orbits. 



so that the stream density in configuration-space is proportional to 
the inverse of this determinant: 



Pstream 



det (D_(t; xo] 



(23) 



At time to the configuration-space distortion tensor equals unity. 
Thus, if we norm the stream density to its initial value, we get the 
following relation for the normed stream density: 

1 



normed 

Pstream V / 



det (D_{t] Xo] 



(24) 



In the rest of this paper we will almost always discuss this normed 
stream density. 

In Fig. |2] we sketch the distortion of the infinitesimal cloud 
around the reference DM particle. Note the difference between the 
stream density evolution in configuration-space and in phase-space. 
The volume of the small cloud grows in Fig. |2] and is not constant 
anymore! This is a result of the projection from phase-space to 
configuration-space. Nearby trajectories spatially diverge in time. 

This formalism can be used to identify caustics in a very ef- 
ficient way. If the DM particle passes through a caustic det(_D) 
passes through zero, and the stream density goes to infinity (for per- 
fectly cold DM). A small cloud surrounding the reference particle 
turns inside out as it passes through the caustic. This corresponds 
to a change in sign of the determinant, and thus can easily identi- 
fied numerically. This property allows the location of caustics to be 
mapped accurately even in complex configurations. Note that the 
possibility of sign changes explains why we took the modulus of 
the determinant in Eq. l |23t . 

The complex fine-grained phase-space structure and espe- 
cially the caustics expected in CDM halos are likely to sub- 
stantially enhance annihilation radiation. These effects have 
so f ar been analysed only under quite simplified conditions 
jHogan 2001: Bergstrom et al. 2001;_Pieri & Branchini 2005,; 
Mohavaee & ShandarinI l2006l : iMohavaee et alj l2007l : iNataraianl 



2007h . As a result, it is unclear how strong such enhancement ef- 



fects will be in the proper cosmological context. Previous studies 
of annihilation radiation from N-body halos had realistic forma- 
tion histories but were unable to resolve caustics. They estimated 
emissivities from the local mean CDM densi ty, thus effectively 
excluding contributions from single streams jStoehr et al.l 120031 : 
iDiemand etn]|2007h . iHoganI ( l200ll) noted that this results in an 
underestimation of the annihilation rate, and suggested that anni- 



hilation might, in fact, be dominated by contributions from the ne- 
glected caustics (which he baptised as "micropancakes"). 

Our formalism enables a robust and accurate calculation of the 
contribution to the annihilation radiation from individual streams. 
The annihilation rate for each particle due to encounters with other 
particles in its own stream can be evaluated directly from the lo- 
cal stream velocity distribution and density. Integrating these rates 
along the trajectories of all particles produces a Monte Carlo esti- 
mate of the intrastream annihilation rate for the system as a whole. 
This automatically includes the contributions from all caustics and 
is exactly the fine-scale contribution which is missing from stan- 
dard N-body-based estimates of annihilation luminosities. 

We now discuss briefly how this is done. Given the very small 
primordial velocity dispersion ao of CDM, and approximating the 
initial configuration-space density as a constant po we can write the 
phase-space density around a reference particle at the initial time to 
as follows :0 

f {x,to) = poNoexp{-l-{x -X {to;xo)yWo{x -X {to;xo))), 



where 



Wo = a(r'diag(0, 1), 



and TVo = {2n)^^^'^aQ^ . Note that this represents a Gaussian dis- 
tribution in velocity-space and a constant density in configuration- 
space. 

Using S{t) = D(t)5o we obtain the phase-space density 
around the particle at the later time t: 

f{x,t) = poNoexp{-^{x-x{t;xo)yW{t){x''x{t;xo))), 
where 

The configuration-space density around the reference particle at 
x{t;xo) is simply the integral of the phase-space density over all 
velocities evaluated at x[t; xo): 



P(t) = po -3 , 



(25) 



where the velocity dispersions ai{t) are given by l/^\i{t) and 
\i{t) are the eigenvalues of the velocity submatrix of W(t). The 
velocity distribution in the principal axis frame of the velocity el- 
lipsoid centred on the particle's position and velocity is given by: 

=iV(t)exp(-Vdiag(ai(t), cra^, a-^t)^ v), 

where N{t) = l/{{2Tvf^'^ai{t)a2{t)a3{t)). Note that this veloc- 
ity distribution is normalised, i.e. its integral over velocity space is 
unity. 

These quantities allow us to calculate the instantaneous anni- 
hilation rate at each point on the particle's trajectory: 

^^PMf A\oj,{v)vg{v) = P^^^^ (26) 
at m J m 

where m is the particle mass and crA(f) the annihilation cross- 
section. We note that in many WIMP models the annihilation cross- 
section is inversely proportional to encounter velocity, and in this 
case {(Jav) is independent of the actual local velocity distribution 



t denotes the transpose of a matrix 
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jjungmanetaLll 19961) . An image of the system in annihilation ra- 
diation can be constructed by integrating all particles forward over 
a short time interval and summing their annihilation contributions 
into a pixel array. We note that equation ( I26t exhibits near-singular 
behaviour as particles pass through caustics and as a result special 
care is needed to obtain the correct contribution to the annihilation 
luminosity in this situation. This will be discussed more fully in 
later papers. 



3 THE DAMAFLOW CODE 

We have developed the code DaMaFlow to test our GDE scheme 
by analysing the behaviour of streams in a broad range of static 
potentials. DaMaFlow integrates the equations of motion and in 
parallel the GDE for a single orbit within user-specified potentials. 
This requires solving 3 + 18 second-order differential equations 
in parallel. The integration algorithm can be chosen to be a sym- 
plectic second-order leapfrog (Drift-Kick-Drift or Kick-Drift-Kick 
formulation) or alt ernatively the ene rgy-conserving and adaptive 
Dop853 algorithm dHaireretalJI 1993 1 of order 8 that allows dense 
output and is very fast. Studies focusing on complex phase-space 
structures, especially in the field of chaos analysis, often use the 
Dop853 algorithm (or even higher order schemes) because of its 
high precision. On the other hand N-body codes often implement 
the leapfrog method because it is the best compromise between per- 
formance and accuracy. We find that with a moderate fixed time- 
step, both formulations of leapfrog are able to give comparable re- 
sults to Dop853. This is an important point because it is not possible 
to run N-body simulations with slow but accurate high-order ODE 
solvers like Dop853. 

DaMaFlow is also set up to do a Numerical Analy- 
sis of Fundamental Frequencies (NAFF) jLaskaJ Il988l Il990l : 
IPapaphilippou & Laskar" l996l : lLaskaij|200 3') and of the resonances 
associated with the chosen orbit. The fundamental frequencies are 
revealed by an integer programming routine. This is needed so that 
we can study the relation between the existence of well-defined fun- 
damental frequencies or resonances and stream density behaviour. 
The NAFF method determines a quasi-periodic approximation to 
the orbital motion. For ordinary Fast-Fourier-Transforms (EFT) the 
accuracy of the determination of the frequencies is of the order of 
1/T, where T is the sampling interval. The NAFF method has an 
accuracy of l/T*. Thus it makes spectral analysis a lot faster com- 
pared to classical Fourier techniques, for example, those used in 
iBinnev & Spergell ( fT982i) . 

To scan large parts of phase-space, DaMaFlow can be run in 
parallel on batch systems in order to integrate a large number of 
different orbits simultaneously. A fast automated stream density fit- 
ting routine was built in to facilitate efficient analysis of the under- 
lying phase-space without user interaction. Before this fitting can 
be done, the stream density has to be smoothed to remove the large 
density spikes produced by caustics. DaMaFlow does this by ex- 
tracting and fitting the lower envelope of the stream density time 
series. This envelope is constructed while the orbit integration is 
running by an iterative on-the-fly minimum finder. 

From a numerical point of view, solving the GDE is quite dif- 
ficult in chaotic regions of phase-space. In these regions the in- 
finitesimal phase-space volume around the reference particle gets 
distorted very strongly. This produces large numerical values in the 
phase-space distortion tensor. And this can lead to overflows and 
round-off errors in numerical computations. In chaos analysis it is 
an established method to do some kind of re-norming to suppress 



these problems. We can do something similar to follow the evolu- 
tion of phase-space density evolution. We use the transitivity of the 
phase-space distortion tensor: 

Dti^ts = Dt2^t3 Dti-.t2 (27) 

where Dt--,tj with i < j is the solution of the GDE with 

Dti~,t (ti) — 1 evaluated at time tj. So the phase-space volume 
can be written as: 

det (^Dt.^t:,^ = det (pt^^t:,^ det (Ot^-^t^) (28) 

Thus dividing the time integration interval and re-initialising the 
distortion after each interval avoids large numerical values. This 
is very similar to the re-norming techniques used for calculating 
the largest Lyapunov exponents of chaotic systems, where the re- 
norming frequ ency is chosen to be o f the order of the dynami- 
cal time-scale jLichtenberg & Lieberm an 1983; El-Zant 2002). Al- 
though this approach works nicely for such applications, it does 
not help us when calculating the stream density evolution, because 
we need the entire phase-space distortion information from ini- 
tial to final time. It is not possible to separate the configuration- 
space part of the phase-space distortion and to do a similar re- 
norming. Thus one cannot avoid large numbers during the cal- 
culation. DaMaFlow therefore calculates all quantities in double- 
precision (64 bits=8 bytes). Even in chaotic regions this is enough 
to follow the system for a substantial amount of time. We note that 
the phase-space density calculation involves the determinant of a 
6x6 matrix, whereas the stream density only involves the deter- 
minant of a 3 X 3 matrix. As a result stream density calculations 
are less strongly affected by large numbers and overflows. We note 
that special software libraries can provide even higher precision 
(e.g. the GMP libraryQ). 

Since we wish to implement the GDE formalism also into an 
N-body code, execution speed and memory consumption are im- 
portant considerations. For each particle we need to store the full 
6D phase-space distortion tensor and the particle's tidal tensor. This 
results in 36 + 6 extra numbers per particle. Thus a 500"^ particle 
simulation in double-precision needs already about 39 Gbytes of 
random access memory (RAM) just for the GDE calculation, as- 
suming we store all information for every particle. 



4 INTEGRABLE POTENTIALS 

Before presenting some results for the evolution of stream densities 
in integrable potentials we want briefly to discuss the choice of the 
initial sheet orientation that goes into y_ (xo) in equation l |19l l. This 
depends, of course, on the starting time to and the problem that is to 
be studied. For example, in a cosmological context V_ (xo) is given 
by the linear initial conditions, so by the coupled initial density and 
velocity fields. We note that, to zeroth order, the sheet orientation 
is simply given by the Hubble Flow: v[xq, to) = H{to)Xj^. 

What is the impact of the initial orientation on later evolu- 
tion? In Fig.[3]we show the evolution of the normed stream density, 
as calculated from eq. l l24b . for a single orbit with four different 
choices of initial stream orientation. For this test we have used a 
spherical Hemquist potential: 



^ http://www.gmplib.org 



The fine-grained phase-space structure of Cold Dark Matter halos 7 




Figure 3. Stream density evolution for different initial sheet orientations 
V_ (xq) in a static Hemquist potential. The general shape of the four curves 
is very similar. The long-term behaviour does not in general depend on the 
initial sheet orientation. 



with M = 1.86 X IO^^Mq and a = 34.5 kpc. The reference 
particle begins its orbit with Xj^ = (35, 17.6, 4) kpc and ?;q = 
(—316.9, —16.3, —4) km s~^. The initial sheet orientations were 
chosen to be (in units of km s^^ kpc^^): a : 0, h : 1, 




and d : 



1 

-1 
2 



It is striking that all four curves have very similar shape and 
very similar caustic spacings, although the caustic locations vary. 
The long-term behaviour does not depend on initial sheet orienta- 
tion, at least in this case. Orientation c produces lower densities 
than the others because the scale of V_ (xo) is larger, but the shape 
of t he lower envelope is v ery similar. 

iNataraian & Sikivig ^2006) show that the caustic shape in 
configuration-space is, in general, affected by the relative size of 
the V_ (^o) matrix elements. A detailed analysis of caustic shape 
thus requires choosingj/ (xg). For example, in their model for 
the Milky Way halo Nataraian & Sikivid ( I20061) initialised trajecto- 
ries at the turnaround sphere with a V_ (aJo) loosely motivated by 
tidal torque theory. This restricted the form and scale of the matrix, 
but still left a lot of freedom. Here our main motivation is not to 
analyse caustic shapes, but rather the long-term behaviour of the 
fine-grained phase-space distribution, in particular of stream densi- 
ties. The initial sheet orientation is thus not an important issue for 
us. In the following we will consider quite general orbits, but will 
arbitrarily set Y_ (xo) — unless otherwise stated[f|. Note that the 
choice of V_ (xo) does not influence the dynamical evolution of 
the distortion tensor as it is not part of the initial conditions for the 
GDE. Only the final projection to configuration-space is affected 
by initial sheet orientation. 

From a dynamical point of view, static, integrable potentials 
are very simple. The motion within them can be described in terms 
of action/angle variables and their Hamiltonian can be expressed 
solely as a function of the actions H = H{J). All motion in these 



potentials is regular, so there are no chaotic regions in their phase- 
space. In action-angle space the orbits lie on tori and are charac- 
terised by a fixed number of fundamental frequencies. DM particles 
in integrable potentials will experience on ly phase mixing. 

Because of these simple properties iHelmi & White! ( 1 19991) . 
hereafter HW, were able to develop an analytic linearised treat- 
ment based on action-angle variables to derive results for the stream 
density behaviour. In their paper they did not specifically focus 
on CDM, but rather analysed how Gaussian clouds in action-angle 
space evolve with time. 

As an example of an integrable potential, we apply our method 
to several Eddingto n potentials ^{r,9) — $i(r) + r){fi cos 9) / 
( lLvnden-Belllll962h . These are separable in spherical coordinates. 



The third integral for this type of potential is /a 



+ 



ri{f3 COS 9). We chose the following example of an Eddington po- 
tential, 



$(r,e) 



yn2 2 

2 , / 2 , ,2^ , P COS 

log (r + d ) H ^ 



(30) 



with Vh ~ 1 km s ^, d = 1 kpc, 13 = 2.5 kpc km s ^ Q and 
studied an orbit which starts at Xjj — (5, 3, 2) kpc with a velocity 
of Uq = (0.62,0.62,0.104) kms~\ 

In Fig.|4]we show the evolution of the stream density for this 
orbit, i.e. the projection from phase-space to configuration-space 
for an initial condition with V_ (xq) = 0. Here and elsewhere (un- 
less otherwise stated) we define 'the 'orbital period" as the radial 
oscillation period for the purpose of making such plots. The late- 
time behaviour of the stream density can be fitted by an analytic 
formula derived by HW: 



p{x, t) = A 



sin 9\prPe\ (t/*orbitai)^ 



(31) 



with only one fitting parameter A. Comparing this to the result for 
the long-term behaviour in HW (Eq. 37) it is clear that A just re- 
flects the initial phase-space distribution and the orbital parameters 
(the derivatives of the fundamental frequencies with respect to the 
actions). The results in HW are calculated for a Gaussian cloud in 
phase-space, not a cold sheet as in our case. Since the initial distri- 
bution only affects A, the long-term behaviour of the two configu- 
rations is the same, as shown in Fig.|4] 

One can clearly see that our method produces caustics at the 
correct positions and is also able to recover the secular evolu- 
tion of the stream density, the l/(t/torbitai)^ density decrease. We 
note also the initial quasi-exponential stream density decrease that 
is often referred t o as M iller's instability (Miller 1964). Recently 
iHelmi & Gome3 feOOTh . hereafter HG, showed this is a generic 
feature of Hamiltonian dynamics an d not, as long believed, an ar- 
tifact specific to N-body integrations jHemsendorf & Menittl200^ : 
iKandruD & Siderisl2003h . HG discuss the effect in detail for spher- 
ical potentials. All our tests with spherical, axisymmetric and triax- 
ial potentials show a similar quasi-exponential initial decay. Since 
DaMaFlow integrates the equations of motion for a single particle 
in a perfectly smooth potential, it is evident that this behaviour has 
nothing to do with N-body effects. 

A NAFF frequency analysis of particle orbits in the Eddington 
potential reveals, as expected, three linearly independent frequen- 
cies. It is this number that dictates the speed with which stream 
densities decrease in static, separable potentials. One can see this 



^ Orbits stalling on axes of symmetry in phase-space can show non-generic 

stream density behaviour for y_ (xq ) = as we will discuss in the section ^ These values do not have any specific meaning. We have chosen them 

on non-integrable potentials (section 5). just in a convenient way. 
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Figure 4. The stream density evolution calculated using the GDE integrator 
DaMaFlow is compared to the a nalytic result obtained b y a linearised treat- 
ment in action/angle variables i Helmi &WhitdlT99i) . The results agree 
essentially perfectly. Notice how well the numerical calculation matches 
the caustics. The upper panel clearly shows that the numerical result also 
has the correct l/(Vorbital)^ behaviour at late times. The initial quasi- 
exponential decay is also visible. 

very clearly from Fig. |5] Here we show the density decrease in a 
simple Kepler-like toy-model: 



$(r) = . 



(32) 



For a — \ the orbit was started from Xq = (—0.33,0.97,0) 
and «g — (—1.0, —0.07, 0) with an energy E — —0.5. This or- 
bit has only one fundamental frequency. Changing a to 0.75 in- 
creases the number of frequencies to two; the loops of the orbit no 
longer close. The starting point for this second case has been set to 
Xg = (-0.27, 1.28,0) andug = (-0.76,0.24,0), corresponding 
also to E = —0.5. The increased number of frequencies results 
in a more rapid decrease in stream density: l/(f/torbitai) for one 
fundamental frequency, and l/(t/torbitai)^ for two. In this sense 
the long-term behaviour of streams in static, integrable systems is 
very simple and is determined only by the number of fundamental 
frequencies. 



5 NON-INTEGRABLE POTENTIALS 

Analytic methods like the HW formalism are only able to deal with 
integrable potentials. This is a serious limitation since realistic po- 
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Figure 5. Stream density evolution for the normal Kepler potential (a = 1) 
and for a modified potential (a = 0.75). Orbits in the Kepler potential 
have a single fundamental frequency. Changing the potential exponent from 
a = 1 to a = 0.75 increases the number of fundamental frequencies 
to two. As a result, the long-term stream density behaviour changes from 
l/(Vorbital) to 1/ (< /<orbital )^ ■ Fof integrable potentials the long-term 
stream density decrease on each orbit is dictated by the number of funda- 
mental frequencies. 



tentials are rarely integrable. To demonstrate that the GDE method 
can also deal with more complex phase-space structure we now dis- 
cuss the well-known ellipsoidal logarithmic potential with a core. 



1 2 

'^{x,y,z) = -vo In 



I 2 , 2 



+ {y/qf + {z/pf) , (33) 



analysing its stream density behaviour and its phase-space struc- 
ture. There are two reasons why we have chosen this potential: 
first there has been substantial previous work on its phase-spac e 
structure teinnev&Spergelll 19821 : IPaDaphiUppou & Laskailll998h . 
so we can compare directly with these earlier results. Second this 
potential is often considered as a good model for galactic halos 
because it reproduces a flat rotation curve and its shape can eas- 
ily be tuned by two parameters that correspond to the axial ratios 
of the ellipsoidal isopotential surfaces. It has been used, for exam- 
ple, for dynamical studies of the debris streams of the Sagittarius 
dwarf galaxy (Helmi 2004). The ellipsoidal logarithmic potential 
without a core (r-c = 0) has also been used to study the influence of 
halo shape on the a nnual modulation signal in dark matter detectors 
tevans et alj200(]h . 

It is known that, depending on the degree of triaxiality, the 
phase-space of the logarithmi c potential can be occupied to a 
large extent by chaotic orbits dPapaph ilippou & Laskailll998h . In 
Fig. [6] we show how the stream density evolves along one of 
these chaotic orbits. This orbit was integrated in a potential with 
q — 1.5, p = 0.5, vo = l,rc = 1,E — 3 and started at 
Xg = (10, 1, 5), = (0.16, 80, —0.16). Here we apply the same 
system of units as IPapaphilippou & Laskad ( 1 1 9981) and write all 
quantities as dimensionless numbers. It is evident from this plot 
that the system mixes very rapidly along this orbit. This is chaotic 
mixing, and contrasts with the phase mixing that we found before 
for regular motion in separable potentials. We note that chaotic or- 
bits are difficult to handle from a numerical point of view because 
of the rapid spread in phase-space that characterises them. To check 
whether we can rely on our stream density values, we also calcu- 
lated the 6-D phase-space density along this orbit. Over the full 
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Figure 6. Stream density evolution along a chaotic orbit in the ellipsoidal 
logarithmic potential with q = 1.5, p = 0.5. The density decreases very 
rapidly, reflecting the chaotic mixing along this orbit. Note that the decrease 
is no longer a power-law, as in the case of regular motion. 

integration range shown in Fig. [6] (40 orbits) it remained constant 
to an accuracy of 10^^. 

Using DaMaFlow we have scanned the phase-space of box 
orbits within the logarithmic potential, identifying chaotic and reg- 
ular regions by calculating the stream density evolution for about 
5 X 10* different orbits. Each orbit was integrated for a fixed time 
interval of 2000 using 2x10''' time-steps. This corresponds to about 
10^ orbital periods. We chose this very long integration time in or- 
der to distinguish between chaotic and regular regionstf]. 

Our results can be compared directly to previous work 
where the same potential was analysed using frequency shifts 
jPap aphilippou & Laskar 1998). This method is based on the fact 
that chaotic motion, contrary to regular motion, has no stable fun- 
damental frequencies, so that frequency estimates shift with time. 
By looking for such shifts one can distinguish between chaotic 
and regular motion. Fig. |7] shows maps of orbit type for two dif- 
ferent sets of axial ratios q,p. This figure can be compared di - 
rectly to Fig. 6(c) and Fig. 4(b) in lPapaphilippou & Laskaj jl998l) . 
For the first of these calculations we adopted the following values: 
q = l.8,p = 0.9, vo = V2,rc = 0.1 and E = -0.404858. For 
the second case we changed the potential shape by instead taking 
g = 1.1, p = 0.9. 

The maps of Fig. [7] are constructed as follows. We start each 
orbit at the centre of the potential to be sure to get a box-orbit with 
zero angular momentum. Each individual orbit can be labelled by 
its initial and Vy velocity components. The value of is then 
determined by the chosen value of the energy, E — —0.404858. 
This is, of course, a very special point within the potential, and 
it turned out that our standard choice of initial stream orienta- 
tion, V_ (xo) = 0, produces highly non-generic stream density 
behaviour in this case; the stream density remains constant! In 
order to get properly representative behaviour we therefore took 
V_ (xo) = 1 when producing Fig.|7] With this set-up, we scanned 
the whole — Vy plane and saved the stream density of each orbit 
after a fixed amount of time. The greyscale in the plots corresponds 

* For chaotic orbits elements of the distortion matiices can become very 
large. Here we are only interested in separating chaotic and regular orbits 
reliably. Thus we do not care about round-off errors in the calculation of 
these matrices. 



to the stream density decrease after that fixed time. Black points de- 
note a very rapid stream density decrease, thus regions of chaotic 
mixing. Grey regions show a slower density decrease, reflecting 
phase mixing and regular motion. Much of the box-orbit phase- 
space is chaotic for q — 1.8, p — 0.9. Reducing the asphericity 
to q = 1.1, p — 0.9 results in a much larger fraction of the boxes 
being regular. 

A comparison of Fig. [7] to Fig. 6(c) and Fig. 4(b) in 
IPapaphilipDou & Laskaj ( Il99& ) shows excellent and detailed 
agreement. Regions of high frequency shift correspond, as ex- 
pected, to those of rapid stream density decrease and chaotic 
mixing. Thus the GDE and the frequency shift technique work 
equally well for delineating regions of chaotic and normal 
phase mixing. We note that the Lyapunov exponent technique 
for identifying chaotic behaviour is closely re lated to the GDE 
jLichtenberg & Liebennarilll983l : lErZantll2002l) . since these expo- 
nents are obtained from the eigenvalues of the 6D distortion tensor. 
Identifying and characterising chaotic behaviour is not the main 
goal of our work here, so we will not pursue this connection further 
in this paper. 

So far we have classified orbits as either regular or chaotic, 
but the regular part of phase-space has substructure in the form of 
resonances ( Carpintero & Aguila3 [T998l : IWachlin & Ferraz-Melld 
iMerritt & Valluriijil9 99). These are regions where the fre- 
quencies of motion are commensurate mioji + m2UJ2 + mscj^ 
where the rrii are integers and the LJi are the three frequencies 
of the regular motion. As shown above, resonan ce influence the 
stream density behaviour over long timescales ( iHelmi & Whit3 
ll999l : ISiegal-Gaskins & Valluril 20071) . This is because they restrict 
the motion to a lower dimensional region in phase-space, implying 
that the orbit does not fill its KAM (Kolmogorov-Arnold-Moser) 
torus densely. In simple terms, the system cannot spread as fast as 
it would do in the non-resonant case because it occupies a space of 
lower dimension. 

Fig. [8] shows the stream density evolution along three dif- 
ferent box orbits for q = 1.8, p = 0.9, wo = V2,rc — 
0.1, E = -0.404858. The initial conditions for these orbits 
were chosen so that they have different numbers of orbital res- 
onances (non-resonant, one resonance, two resonances): = 
(-0.08, -0.70, -0.090), = (-0.60,-1.50,-0.07) , x^ = 
(-0.01,-0.67,-0.08), = (-0.53,-1.60,-0.06) and x^ = 
(0.08, -0.63, -0.14), = (-0.51, -1.51, -0.52). It is clear 
that the number of resonances has a major effect on the final stream 
density decrease over timescales similar to those shown in this plot. 
The difference in stream density between the non-resonant case and 
the case with two resonances is about 3 orders of magnitude after 
350 orbits! Resonances also have a strong influence on the shape 
of the orbit in configuration-space, as shown in the lower panel of 
Fig. [8] Two resonances restrict the orbit to a line. From the shape 
of the orbits it is evident why the stream density changes so much 
with the resonances. The particles cannot spread over a large region 
if they are confined to a space of small dimension. 

The regular region of phase-space for the logarithmic poten- 
tial is occupied by resonant and non-resonant regions. We used 
DaMaFlow to scan the regular region (as for the chaos maps above) 
and fitted the stream density decrease by a power law in time. 
A non-resonant motion then gives an exponent of 3, whereas re- 
gions with two resonances should give 1. We binned these power 
law exponents (bin size 0.1) and coloured them. The result of this 
procedure is shown in Fig. |9] For this map we integrated a total 
450 X 450 orbits for about 10'* orbital periods. Each integration re- 
quired 25 X 10*' KDK leapfrog time-steps. The chaotic regions are 
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Figure 7. Chaos maps for box orbits in the logaiithmic ellipsoidal potential 
for two different sets of axial ratios. The upper plot is for a highly aspherical 
potential with q = 1.8, p = 0.9, whereas the lower plot is for a rounder 
potential with q = 1.1, p = 0.9. The greyscale indicates the stream den- 
sity decrease after a fixed time interval (2000 time units or about 100 orbital 
periods). Black regions mark the very low final stream densities found for 
chaotic orbits, while grey regions mark the higher stream densities found for 
regular orbits. Densities decay quasi-exponentially in the former case, but 
only as a power law of time in the latter. These plots can be directly com- 
pared to Fig. 6(c) and Fig. 4(b) in Papaphilippou & Laskat ( 1998), where a 
frequency shift analysis of the same system reveals exactly the same struc- 
tures. 



shown in grey, while the regular regions are shown in three different 
colours depending on their stream density behaviour. We note that 
the chaos detection here was carried out by imposing a threshold 
10~^® on the stream density decrease. Every orbit that crosses this 
threshold during its evolution is considered as chaotic and marked 
as grey in the map. This is the reason why the chaotic pattern here 
is not identical to that in Fig.|7] 

Most of the regular phase-space in this figure is occupied 
by non-resonant orbits shown in green. Superposed on these is 
a fine network of resonance lines shown in red. In blue regions 
the stream density decreases linearly with time because there 
are two resonances. Note how well our method locates the res- 
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Figure 8. Stream density evolution on resonant orbits. Stream density drops 
at a rate which depends on the number of independent orbital frequencies. 
Non-resonant orbits have three independent frequencies, and their stream 
density decreases like (i/*orbital)~^ at late times. Resonances reduce the 
number of independent frequencies and slow the decrease of stream densi- 
ties. With one resonance there are two independent frequencies; the stream 
density then falls as the inverse square of time. With two resonances, only 
one independent frequency remains and density falls as the inverse first 
power of time. The number of resonances also strongly affects the orbit 
shape in configuration-space. As is visible for the examples in the lower 
plot, non-resonant orbits (red) fill a 3D volume, whereas orbits with two 
resonances (green) are restricted to a line. 



onance li nes in the initial conditi on space spanned by Vy and 
Vx- Papap hjrippou & Laskad ( 1 19981) analysed resonances with fre- 
quency maps by plotting the rotation numbers defined as ai = 
vl/vs and 02 = vm/vs, where Vi are the fundamental frequen- 
cies along the long (L), short (S) and middle (M) axes. It is straight- 
forward to identify the resonance lines in Fig. |9] with those in the 
frequency map of Papaphilippou & Laskar (1998) by applying a 
NAFF frequency analysis to the orbit corresponding to any spe- 
cific initial condition, for example, one on a given resonance line. It 
turns out that we can identify all resonance lines in Fig.|9]with sim- 
ilar lines in the frequency map. At the intersection of these lines we 
have periodic orbits (satisfying two reson ance conditions) with the 
same rotation numbers as those found in IPapaphilippou & LaskaJ 
l ll998l) . For example, the line going from the upper left corner to 
the lower right corresponds to the (3, 1, —2) resonance, meaning 
that 3ai + la2 - 2 = 0. 
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Figure 9. Resonance structure of a logarithmic potential as revealed by 
the GDE. Green indicates regions with no resonance (i.e. three indepen- 
dent orbital frequencies). Stream densities for these orbits decrease as 
(Vorbital)"'^ at late times. Red indicates orbits with one resonance, 
for which stream densities decreases as (t/torbital)^^- Orbits with two 
resonances are coloured blue. This map was constructed by integrating 
450 X 450 orbits, each for about 10* orbital periods, using DaMaFlow 
in parallel. Each orbit required 25 X 10*' time-steps using a KDK leapfrog 
algorithm. 



We conclude that our method can resolve the structure of 
phase-space equally as well as the standard frequency mapping 
technique. 



6 TRIAXIAL DARK MATTER HALOS 

In CDM cosmologies dark matter halos are not spherical. Further- 
more, simulations suggest that their shape should vary with radius, 
both equidensity and equipotential surfaces being rounder (on aver- 
age) at larger radii. Several studies have tried to constrain the shape 
of the Milky Way's halo by analysing the properties of observed 
tidal s treams like t hat of the Sagittarius dwarf galaxy (Ibata et al. 
1200 ll : lHelmill2004t) . Recently, iHavashi et alj j2007l) analysed the 
radial variation in potential shape of simulated halos that might 
correspond to that of the Milky Way. Although there is substantial 
object-to-object scatter, on average they found a relatively rapid 
transition from aspherical to almost spherical which occurs near 
the scale radius Ts of the best fitting NFW profile. They provide a 
simple fitting formula for this mean behaviour. 



log ( — or — ) — a 

a a 



tanh ( 7 log — 



1 



(34) 



for the principal axial ratios b/a and c/a. (Note that they actu- 
ally provide two different sets of fitting parameters for Eq. ([34} 
depending on the principal axial ratios.) They also propose a mod- 
ified NFW potential that takes into account the variation in shape, 
but this potential is not very convenient because it is not straightfor- 
ward to derive the corresponding equations of motion. The exam- 
ples given above show that potential shape can have a substantial 
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Figure 10. Isopotentials for the outer and inner parts of one of our triax- 
ial NFW halos. It is obvious that the halo becomes rounder as one moves 
outwards. In this case the transition scale Ta was chosen to be equal to the 
scale radius of the NFW profile. 



effect on stream density evolution, so it is interesting to see how 
strong such effects can be in a realistic model. 

To analyse this we have built a simple extension of the NFW 
model that qualitatively reproduces the shape variation found by 
Havashi et al. (2007) but which has simple equations of motion 
that can easil y be implemented in DaMaFlow. (For another simi- 
lar model, see lAdams et al.l ( |2007|) .) 

We model the variable shape of the NFW halo by replacing 
the euclidean radius in the formula for the potential of a spherical 
NFW halo by a more general "radius" f given by: 



(rg + r) rs 
{ra + t-e) 



(35) 



Here is a transition scale where the potential shape changes from 
ellipsoidal to near spherical and r_B is an ellipsoidal "radius" given 
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Figure 11. Com parison of the radial variation of isopotential axial ratios 
for N-body lialos jHavashi et alj|2007l) to that for our simple triaxial NFW 
model. The N-body values are the average of those found in lHavashi et alj 
<2007l) . The transition scale Va of our model has been set equal to the scale 
radius Ts of the underlying NFW profile. 



by: 



(36) 



where we require + 6^ + = 3. Thus for r ^ ra f = re and 
for r 3> T'a r = r. We then take the potential to be ^{x, y, z) = 
^ NFw(^(a^;?/i ■?)) whic h reproduces the general behaviour found 



bv lHavashi et alj 1 2007h with a smooth transition around ; 



For a specific example, we have chosen the transition scale to 
be the scale radius of the NFW profile and have taken values for 
a, b and c that give central principal axial ratios that are compara- 
ble to those found by Hayashi et al: b/a = 0.78 and c/a = 0.72. 
Our choice is a = 1.18, b = 0.92, c = 0.85. For the NFW profile 
we used a concentration of r2oo/rs = 7.0. We checked Poisson's 
equation for this potential to ensure that it implies a positive den- 
sity everywhere. The check was performed by DaMaFlow evalu- 
ating the negative of the trace of the tidal field on a fine 3D grid. 
This is just the Laplacian of the potential and so proportional to the 
corresponding density. Since the density field is continuous, posi- 
tive density values on the grid should guarantee a positive density 
everywhere. 

Fig. [To] shows isopotentials in the outer and inner parts of 
the halo. All distances are expressed in terms of the scale radius 
Vs of the NFW profile. The transition from spherical to aspherical 
can clearly be seen as the centre is approached. Fig. 1 11 [compares 
the radi al variation in axial r atios in our model and in the simula- 
tions of iHavashi et alj j2007h . The simulation axial ratios are cal- 
c ulated with eg. 34t u sing the average values for a, 7, found 
in lHavashi et al.l 1 20071) . The lines for our model are calculated as 
follows. For a given value of f we computed the intersections of 
the corresponding isocontour with the x, y and z axes. So we get 
three values a" , b" , c" . To look for their variation over distance we 
define the mean distance r" = \/ (a ")^ + jb")'^ + (c")^ . This is 
essentially the same procedure which lHavashi et al] ( l2007h applied 
when fitting the isopotentials of their simulated halos. Thus we can 
compare directly with their results as in Fig.[TT] The qualitative be- 
haviour of our model is very similar to that of the simulations. It 
is not necessary to dema nd an exact fit since the scatter between 
different halos studied bv lHavashi eTZI ( l2007h is quite large. 

We implemented this potential into DaMaFlow and looked at 



four different orbits with the following apo-/pericentre distances 
in units of r^: 8.9/6.1, 20/5.9, 2.2/0.5, 1.2/0.4. We compared the 
stream densities predicted for our triaxial model to those predicted 
in the corresponding spherical NFW profile. We fixed the starting 
point and the velocity direction to be the same in the two cases. 
The amplitude of the velocity was then set to give the same energy 
in the two cases. With this procedure the orbits covered compa- 
rable regions in configuration-space and had nearly the same apo- 
and pericentre distances. Note that it is impossible to get identical 
orbital shapes in the two potentials. Especially in the inner parts 
of the halo, where the two potentials differ substantially, the orbits 
have different shapes. In a spherical potential orbits are confined to 
a plane by conservation of angular momentum, but this is not the 
case in a triaxial potential. 

In Fig.[T2]we show the stream density evolution for these four 
orbits. Two belong to the outer halo with peri- and apocentre be- 
yond the scale radius. As expected their stream density behaviour 
is very similar in the two potentials. As soon as orbits penetrate 
the inner halo, however, the behaviour is quite different in the two 
cases. After 100 orbits, streams are about 100 times less dense in 
a triaxial halo than in a spherical one. Note that this is just what 
one would predict, given the (t/torbitai)~^ and (f/torbitai) den- 
sity evolution expected for regular motion in spherical and triaxial 
potentials, respectively. 

In the case of a cored ellipsoidal potential we showed above 
that, depending on the level of triaxiality, substantial fractions of 
phase-space can be occupied by chaotic orbits. Thus we may ex- 
pect such orbits to be present in our triaxial NFW profile also. 
Previous work on galactic dynamics has demonstrated the pres- 
ence of chaotic orbits in the potential s corresponding to a va- 
riety of cus py, triaxial density profiles ( Merritt & Fridman 199] 
Merritt & Valluri 1996 : Valluri & Mer ritt 1 998 ; Ka ndrup & Siotjiu 
200 3I : Capuzzo-Dolcetta et alj|2007l) . To search for chaotic orbits 
in our model, we have integrated 2 x 10* different representative 
orbits and studied the predicted behaviour for the density of their 
associated streams. For simplicity we have chosen the initial con- 
ditions for these orbits at random from the known analytic distri- 
bution function of a Hemquist sphere matched to the mean radial 
density profile of our triaxial NFW model. A self-consistent Hem- 
quist sphere has density profile and potential: 
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We matc h to our NFW model by an appropriate choice of the scale 
length a dSpringel et alJl2005h . A standard inversion technique can 
then be used to select a random set of initial orbital positions from 
this (spherically symmetric) distribution. Appropriate initial veloc- 
ities c an be generated by applying the von Neumann rejection tech- 
nique ( jPress et al.|[l992l : lAscasibar & Binnevll2005l) to the analyti- 
cally k nown distributio n function of the self-consistent Hernquist 
model ( iHernquis^l 1 990l) . This procedure does not, of course, sam- 
ple orbits with a weighting which would self-consistently repro- 
duce our triaxial NFW model. Nevertheless, the similarity of the 
NFW and Hernquist models should ensure that our selected orbits 
cover the regions of phase-space which would be populated in a 
truly self-consistent model in a reasonably representative way. This 
is sufficient to evaluate the overall importance of chaotic orbits in 
the model. 

We integrated each orbit as long as required to guarantee 
max(a;cross, J/cross, .Zcross) = 10^, where icross is the number of 
crossings along the i-coordinate, i — x,y, z. The integration was 
done by the Dop853 algorithm with a very high precision to get an 
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Figure 12. Stream densities in a spherically symmetric NFW potential are compared to those expected in a more realistic DM halo with a shape that varies 
with radius. Orbits with pericentre inside the transition scale r'a = Ts show a substantial difference between the two cases. After about 100 orbits streams are 
roughly 100 times less dense in a triaxial halo than in a spherical one. Thus spherical models for the Milky Way's halo are likely to underestimate the number 
of streams in the solar neighbourhood by two orders of magnitude. 



relative energy error below 10~ over the whole integration time. 
We have chosen such high energy conservation to be sure that the 
integration works correctly even though the potential is cuspy. 

We plot the final stream stream density using a greyscale just 
as in Fig.|7] The axes of the resulting plots in Fig.[T3]are the orbital 
energy in the triaxial NFW potential and the circularity based on the 
spherical Hernquist profile that was used to generate the initial con- 
ditions. Black points correspond to very large decreases in stream 
density, hence to non-regular motion. The upper panel shows all 
orbits whereas the lower panel takes only box orbits into account. 
We distinguished box and tube orbits by looking at the angular mo- 
menta around the symmetry axes. For a box orbit the sign of all 
three momenta changes along the orbit. On the other hand a tube 
orbit has one axis along which the angular momentum has a fixed 
sign. 

First of all it is quite striking from these plots that there are 
orbits that are not regular and show up as black points. A line of 
fixed pericentre (l/Sr^) in Fig. [7] clearly shows that these are or- 
bits that reach the innermost part of the halo, where they feel the 
strong triaxiality and the cusp. This is not surprising because pre- 
vious studies showed that cuspy, triaxial potentials exhibit chaos 



jValluri & Merrittlll998l ; iKandrup & SiopisI l2003b . These studies 
found that box orbits are primarily affected but also some tube or- 
bits. This behaviour can clearly be seen in Fig. 1131 

It is not clear how to distinguish between regular and chaotic 
motion based on the stream density decrease after a given num- 
ber of orbital periods. There is no "gap" in the stream density dis- 
tribution which might separate regular and chaotic motion. Meth- 
ods like the freque ncy shift derived from a NAFF analysis run into 
the same problem jValluri & Merritlll998l) . Lyapunov methods are 
better suited to t his problem, but may also have problems making a 
clear distinction jKandrup & Siderisl2002h . As our objective here is 
not a precise chaos analysis of our triaxial NFW potential, we take a 
fiducial value of 10"'^^ for the stream density which separates reg- 
ular and chaotic motion. This value is based on the stream density 
distribution function and corresponds to a stream density below the 
"regular motion bump". Orbits with a stream density below 10~^^ 
will be considered as chaotic. With this criterion about 35% of the 
orbits with binding energy between 0.91 |-Eo| and 1.00 |i5o| are 
chaotic. This fraction is not par ticular high compared to previous 
studies. IValluri & Merritll ( Il998l) . for example, find fractions up to 
80% depending on the parameters they use for their triaxial density 
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tic Centre is roughly 0.2 to 0.3rs. Figure [121 suggests that mix- 
ing will have reduced stream densities from their values at infall 
by roughly 4 orders of magnitude for typical streams in the solar 
neighbourhood. Since the mean DM density near the Sun is sub- 
stantially greater than typical stream densities at infall, we expect 
at least IC CDM streams to contribute to the density in the solar 
neighbourhood. We note that this is still likely to be a substantial 
underestimate, as it is based on a static, smooth halo model and so 
neglects additional mixing effects which may be important, in par- 
ticular mixing in precursor objects, mixing due to scattering by halo 
substructure, and chaotic mixing. Such effects can only be treated 
properly by applying the GDE to structure formation in its proper 
cosmological context. This requires the use of N-body simulations. 
We note that particles orbiting the innermost part of the halo are 
also affected by the disk potential leading to a non-spherical con- 
tribution to the potential. 
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Figure 13. Qualitative view of the phase-space structure of our triaxial 
NFW potential. These plots give results for 2 X 10* orbits. The greyscale 
represents the stream density decrease after max (xcross , J/cross , ^^cross ) = 
10"^. Black points correspond to orbits with a very strong density decrease, 
thus to non-regular orbits. Such orbits are found primarily in the inner re- 
gions of the potential. Orbits with lower binding energy are mostly regular, 
showing a much smaller stream density decrease. The black points were 
plotted last to avoid over-plotting by the grey ones. Eq corresponds to a 
particle at rest at the centre of the potential. The solid red line represents the 
constant pericentre line (r'peri = l/Sr^). In the lower panel only the box 
orbits are plotted. Most of the non-regular orbits are boxes deep inside the 
potential. 



profile. Recent studies of self- consistent models of cuspy tri axial 
galaxies with dark matter halos dCapuzzo-Dolcetta et alj2007l) also 
find chaotic orbits to play an important role. Although we have car- 
ried out only a qualitative analysis it is clear that chaos plays a role 
in the centre of our model also. Box orbits are mostly affected by 
chaotic mixing because they reach the innermost part of the halo. 
We note that the four orbits shown in Fig. [12] are regular. None of 
them has a pericentre distance below the "critical" distance 1 / Sr^ . 

These results show that stream densities near the Sun are pre- 
dicted to be much lower for a realistic triaxial potential than for the 
corresponding spherical potential. The orbital period near the Sun 
is about 2.5 x 10* yrs and the distance of the Sun from the Galac- 



7 FINE-GRAINED PHASE-SPACE ANALYSIS IN 
N-BODY CODES 

The main motivation of our work is the desire to address is- 
sues of mixing and fine-scale structure in full generality by build- 
ing the GDE into current state-of-the-art N-body codes. We have 
done this for the current version of the GADGET code jSpringell 
l2005h . This is a massively parallel N-body code which has al- 
ready been used to c arry out very large cosmological simulations 
dSpringel et alj[200^ . It c alculates the gravitational forces with an 
efficient TreePM method ISode et al.|[200ol : IXij|l995l : lBaglj|2002h 
and uses a domain decomposition scheme based on space-filling 
(fractal) Peano-Hilbert curves to achieve good work-load balance 
in parallel operation. 

To implement the GDE within GADGET we needed to ex- 
tend various parts of the code. The dynamics of the distortion tensor 
are driven by the ordinary gravitational tidal field. The correspond- 
ing tidal tensors for each DM particle have to be calculated using 
the same Tree- or TreePM-Method as the forces in order to provide 
the correct driving term in the GDE. While the forces are given 
by the first derivative of the potential, the tidal tensor is made up 
of second derivatives. The particles in an N-body simulation can 
be thought of as a Monte-Carlo sampling of the real DM phase- 
space distribution, but the coarseness of this sampling introduces 
unwanted discreteness or "two-body" effects which are likely to 
be more serious for higher derivatives of the potential. Such ef- 
fects are usually mitigated by softening the gravitational potential 
of each particle. GADGET uses a spline softening kernel function 
with compact support: 
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The softened potential of a point mass is then given by (x) = 
{Gm/h)W2{\x\/h) with a softening length h = 2.8e, where e 
is the Plummet equivalent softening length. The potential (and so 
force and tidal field) become Newtonian if \x\ ^ h. From this soft- 
ened potential we can calculate the softened tidal field of a point 
mass: 
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Thus the softened tidal field acting on particle k at position (its 
tidal tensor) is given by: 



Sijinigi 



(38) 



+ (a;i,i - Xfc.i) {xij - Xk,j) mig2 



SijTnigi 



+ {xi,i - Xk,i) {xij - Xk,j) mig2 

- SijUikQi (0) 

The last step highlights a difference between the tidal field calcu- 
lation and the normal force calculation. The tidal field is obtained 
using the same tree-walk as the forces. The latter are calculated by 
evaluating the full simi without excluding particle k. This is 
simply to avoid additional bookkeeping; the particle-particle force 
vanishes when I = k so including the self-term does not affect the 
result. This is not the case for the tidal field, for which one must 
add an extra term to the diagonal tidal tensor elements to remove 
the self-tidal field. This is similar to the self-energy correction that 
is needed when using the tree to evaluate the total potential energy 
of the system. 

For larger simulations it is not efficient to use the tree alone. 
In such cases the TreePM method can be much faster. In this 
scheme the potential is split into short-range and long-range parts 
$ = Specifically, in Fourier space GADGET takes 

= $fe exp {-kVi) (39) 

where Vs defines the spatial scale of the force split and should not 
be confused with the scale radius of the NFW profile. The long- 
range potential is calculated by mesh-based Fourier techniques. In 
Fourier space the tidal field can be calculated by just pulling down 
— {ikjY withj = a;, J/, a. The short-range potential in real space is 
given by: 



3>short {x) = i^i) erfc (^^J , (40) 

and the corresponding short-range part of the tidal field by: 
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where is the softened point mass force, xi^i is defined as the 



smallest distance of any of the periodic images of particle I to the 
coordinate Xi of x, and n = ^Jxf ^ + xf ^ + xf ^. 

The time integration also needs modification in order to in- 
tegrate the GDE in parallel with the equations of motion. For this 
it is desirable to write both the equations of motion and the GDE 
in a time-symmetric way. This fits best into GADGET'S quasi- 
symplectic integration scheme which is a second-order leapfrog. 
For the GDE we need to integrate two differential equations of 
second-order to solve for D, with i = xx, xv. Let W_. denote the 
first time derivative of D^. We can then define the system state vec- 
tor § as, 



S=(x,v,D ,D ,W ,W V 

\ ^=XX ^=XV ^='XX ^^xv J 



(41) 



The equations of motion and the GDE can now be written as one 
equation for S: 



§ {t;xo) = f [S {t;xo)) 



(42) 



The right hand side does not depend on the time derivative of S. 
This allows the use of a time- symmetric leapfrog scheme with the 
following Drift- and Kick-operators: 
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(44) 



where T is the ordinary tidal tensor and F_ the gravitational force. 
Although the time integration now needs to solve 18 additional 
non-trivial coupled second-order differential equations, it turns out 
that the loss in performance is not dramatic, even if we do the cal- 
culation for all DM particles in the simulation box. From a com- 
putational point of view, the strongest impact comes from the extra 
memory that is needed to keep track of the distortion tensor. Ev- 
ery particle needs the tidal tensor (6 numbers due to symmetry) 
and the distortion tensor (36 numbers). Nevertheless, with current 
computer capabilities this is not a major limitation. 

In general the state of an N-body simulation is not stored fre- 
quently enough to catch the caustics that occur along each parti- 
cle's orbit. To avoid missing these we implemented a caustic finder 
that examines every drift operation of the time integration. As de- 
scribed above, sign changes in the determinant of the configuration- 
space distortion tensor _D indicate that a particle has passed though 
a caustic. Whenever this happens the event is logged. We are then 
limited only by the time-step of the simulation and this is normally 
small enough to catch all large-scale caustics. 

For flexibility in testing, we have also implemented the GDE 
formalism in a version of GADGET which allows certain static 
potentials, in particular, NFW, Hemquist and cored ellipsoidal log- 
arithmic potentials, to be included in addition to the self-gravity of 
the particles. 

As a first test of our implementation in GADGET, we have 
compared the behaviour predicted for N-body realisations of a 
static Hernquist sphere to that found for an integration in the cor- 
responding smooth potential. To get a system which resembles the 
Milky Way's halo we take M = 1.86 x IO^^Mq and o = 34.5 
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Figure 15. Stream densities for an A*^ = 10^ splierical Hemquist N-body model of the Galactic halo after 5.0 Gyr of integration (blue hatching from bottom 
left to top right) are compared to those found by integrating from the same initial conditions in the corresponding smooth potential (red hatching from bottom 
right to top left). In each case the sphere was divided into nine spherical shells and the normed stream densities of all particles in each shell at the final time 
were histogrammed with bin width 0.5 in log(density). As expected, the lowest stream densities are reached near the centre of the sphere where the dynamical 
time scale is shortest. In the two outermost spherical shells most of the particles have not yet undergone many caustics. The stream density distribution there 
has a strong bias towards values > 1. The agreement between the N-body live halo and the smooth potential is good. 



kpc in equation l l37t . The N-body realisation was constructed as 
described in section 6. 

In Fig. [14] we compare the evolution of stream density for a 
specific orbit in an A'^ = 10^ live Hernquist halo to that found 
when integrating the same orbit in the corresponding smooth po- 
tential. The Hernquist sphere had the parameters given above and 
the particular orbit chosen here had peri- and apocentre of 25 kpc 
and 33 kpc, respectively, giving a period of about 0.5 Gyr It was 



integrated for about 15 orbits or 7.5 Gyr The N-body softening 
was taken to be e = 1.5 kpc. The 6-D phase-space density re- 
mained constant to better than 1 part in 10** in both integrations, 
but the stream density evolution still differs significantly between 
them, in particular in the timing of the caustics and in the detailed 
behaviour of the lower envelope. This is a consequence of the well 
known divergence between nearby orbits in N-body systems which 
is caused by the cumulative effect of many small perturbations due 
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Figure 14. The stream density evolution found by integrating tfie GDE 
along an orbit in a live N-body realisation of a Hernquist sphere is compared 
to that predicted for the same initial condition in the coiresponding analytic 
potential. The N-body halo was realised with A'^ = 10^ particles using a 
softening length, e = 1.5 kpc. The N-body evolution is very similar in 
shape and caustic frequency to that in the smooth potential. The 6-D phase- 
space density remained constant to an accuracy of 10~* over the full N- 
body integration. 



to discreteness jKandruD & Siderisll2003l) . The GDE is very sensi- 
tive to such noise. The features in the two curves are, nevertheless 
very similar, in particular the number and spacing of caustics and 
the overall shape. 

Fig. [15] shows the normed stream densities after 5 Gyr of in- 
tegration for all particles in a live Hernquist halo with TV — 10''' 
and the parameters assumed above. For this integration we adopted 
e = 2.0 kpc. We divide the particles into 9 radial bins contain- 
ing approximately equal numbers of particles and then histogram 
the stream densities, both for the N-body simulation and for in- 
tegrations from the same initial conditions in a smooth Hernquist 
potential. Typical stream densities in Fig. [15] decrease towards the 
centre of the sphere. This is because shorter dynamical times result 
in enhanced mixing in the inner regions. (Recall that stream den- 
sities decrease as (t/torbitai)^^ in a spherical potential, and so are 
smallest where the orbital periods are shortest.) The two outermost 
radial shells are dominated by particles with long orbital periods 
which, as a result, have stream densities of order unity. The high 
stream density tails of the histograms are due to particles which are 
close to caustic passage. They thus have the universal power-law 
shape A^(> p) oc expect ed near a caustic (see, for example, 
iMohavaee & Shandarid (l2006h). 

The two sets of histograms in Fig. [15] are very similar. Al- 
though stream densities evolve differently along orbits from a given 
initial condition in the N-body and smooth potential cases (see 
Fig. I14t the statistical results for ensembles of initial conditions 
are similar. N-body discreteness effects do not cause substantial 
systematic shifts in the stream density distributions predicted for 
this test problem. A small systematic effect is visible at low stream 
densities. The N-body integration produces more very low-density 
streams than the integration in the corresponding smooth poten- 
tial. This is indeed due to discreteness effects, as evidenced by the 
fact that we find the excess to depend on the N-body softening; 
smaller softenings result in a larger tail of extremely low-density 
streams. On the other hand, too large a softening leads to incor- 
rect representation of the mean force near the centre of the system. 
Thus a trade-off is needed to define the optimal softening. This 



has been m uch discus sed with refere nce to conventional N-body 
simulations iMerrittll9 96: Athanasso ula et al.l200ol : lDehnenl200ll : 
iRodionov & Sotnikovall2005i :lZhan 20()6^ but we note that the situ- 
ation is worse for our current application, since the evolution of our 
extended state vector (equationl42b depends on the tidal tensor. The 
additional spatial derivative relative to the force makes our GDE in- 
tegrations substantially more sensitive to discreteness than a stan- 
dard N-body integration. This suggests that the optimal choice of 
softening will be larger for GDE integrations than for conventional 
N-body integrations. 

Fig.llSlshows that the high-density tails of the stream density 
distribution agree well between the N-body and smooth potential 
integrations. This suggests that the number and the strength of the 
caustics must be similar in the two cases. We can check this explic- 
itly by again dividing the sphere into radial shells and then calcu- 
lating the median number of caustic passages by the final time for 
the particles which end up in each shell. In Fig.[T6]we compare the 
results of this exercise for the N-body and smooth potential inte- 
grations using 50 shells. The level of agreement is striking. Only 
within about 3 kpc of the centre is there a significant difference be- 
tween the two curves. This is comparable to the softening used for 
the N-body system, so it is not surprising that particles in this inner 
core pass through fewer caustics in the N-body case. 

The number of caustic passages depends very little N-body 
parameters. In Fig. [17] we plot median caustic count against radius 
for two different mass resolutions and for a fixed softening length 
of 0.5 kpc, four times smaller than in Fig. [T6] After 5 Gyr, the 
highest resolution simulation produces a median caustic count at 
1 kpc which agrees with that for the smooth potential integration in 
Fig. [16] confirming that that the disagreement in that figure was due 
to the softening of the N-body simulation. It is remarkable that par- 
ticle number has no strong impact on the median caustic count. The 
two simulations in Fig.[T7]differ by a factor of 32 in particle mass, 
yet outside 4 kpc they agree very well both with each other and 
with the more softened integration of Fig. [16] The reason for this 
stability is that the caustic count is an integer which is augmented 
only when the determinant of the distortion tensor changes sign. As 
a result, it is much less sensitive to the exact values of the distortion 
tensor elements than is the stream density (which depends on the 
value of the determinant). 

We can estimate the median number of caustic passages for 
particles at radius r very simply as nt/T{r) where t is the age of 
the system, T(r) the period of a circular orbit at radius r, and k a 
proportionality constant. Then r(r) — 27rr/Vc('"), where Vc(r-) = 
\/ GM[r) /r is the circular velocity at radius r. For a Hernquist 
sphere, the mass M{r) within radius r is M{r) = Mr^ j (r + of" . 
Putting all this together we get: 



c(r, t\ k) — K- 



t _ K t^GM/r 



(45) 



T{r) 2n (r + o) 

where c(r, t; k) is the predicted median caustic count at radius r. 
This estimate works very well, with a best fit k — 4.2. The de- 
viation at large radii where the caustic number is low can be ac- 
commodated simply by adding a small constant offset, as indicated 
by the dashed line in the figure. HW (Eq. 37) already showed that 
caustics occur in a spherical potential when pe = or pr ~ 0, 
thus at the turning points in the 9 and r coordinates. If pg and p,. go 
through zero at different times we would expect four caustics per 
orbital period. This is surprisingly close to the value of k that we 
estimate directly from the simulation, given that the particles seen 
at radius r actually have a wide range of orbital periods, rather than 
all having the circular orbit period T(r). 
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Figure 16. Number of caustic passages after 5 Gyr as a function of final 
distance from the centre for 10^ orbits in a Hemquist sphere. The green 
dashed curve gives the median number of caustic passages at each radius 
for pailicles in an N-body realisation of the system integrated with softening 
parameter, t = 2kpc. The red solid curve shows the result when these same 
initial conditions are integrated within the corresponding analytic Hemquist 
potential. The results coincide except within about 1.5 softening lengths of 
the centre. This demonstrates that caustic counting is very robust against 
discreteness effects. 
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Figure 17. Median caustic passage count against radius, as in Fig. 1161 but 
for different particle numbers and for a smaller softening (e = 0.5 kpc). 
The panels show results after 2.5 (top) and 5.0 (bottom) Gyr of evolution. 
Clearly the caustic count goes up between the two times, so the curves are 
higher in the lower panel. The radial dependence is very well represented 
by eq. )45t . Increasing the number of particles by a factor of 32 has very 
little effect on these curves, which also agree with those of Fig. [16] This 
shows that the caustic count along an orbit is very insensitive to the N-body 
parameters (particle number and softening) used to integrate the system. 
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Figure 18. The number of particles that have passed through a given num- 
ber of caustics is plotted against the number of caustics for our A'^ = 
2048000 simulation after 2.5, 5 and 10 Gyr. The thin black lines are an- 
alytic estimates based on eq. (47). 



Rather than focussing on the median count of caustic passages 
as a function of radius, one can examine the distribution of the num- 
ber of caustic passages, i.e. the number of particles that have passed 
through a given number of caustics after some given time. FigurefTS] 
shows such distributions after 2.5, 5 and 10 Gyr for our highest res- 
olution (TV = 2048000) simulation. With increasing time the char- 
acteristic number of caustic passages increases and the number of 
particles with a small number of caustic passages decreases. 

We can make a simple analytic model for these distribu- 
tions based on Eq. l |45t . There are A-k p{r)r^ / m Ar particles in 
the interval (r, r + dr), where m is the mass of a simulation 
particle and p{r) the (analytic) density profile of the Hernquist 
sphere. If we make the approximation that all particles at radius 
r have a caustic count equal to the median count predicted by 
eq. l l45t , the number of particles with caustic counts in the interval 
(c(r-, t\ ft), c(r, t; k) + dc(r, t; n)) will be the same as the number 
of particles in (r, r + dr), so: 



/(c)dc 



1 2 

AiYp(r) — r dr 
m 

_ M r{c,t;K) 1 

47r CL-— 7 -r 

27r (r\c,t; Hi) + a)'^ m 



dr(c, t; k) 



dc 



dc. 



(46) 
(47) 



where r{c,t;K.) is the inverse function of c{r,t;K,). As Fig. 1181 
shows, this formula represents the simulation results very well, sug- 
gesting that the variation in caustic count with radius is more impor- 
tant than the scatter in caustic count at given radius for determining 
the overall shape of the count distribution. 

From these first tests we conclude that our N-body implemen- 
tation is working well, that caustic properties can be predicted very 
robustly, at least when the caustics reflect large-scale structure in 
the system, and that stream densities can also be predicted reliably 
provided care is taken to ensure that discreteness effects are under 
control. 



8 CONCLUSION 

Direct DM detection experiments operate on length-scales far be- 
low the resolution of current structure-formation simulations. The 
fine-grained phase-space structure on these scales will determine 
the signal they see. In addition, small-scale structure can substan- 
tially enhance the annihilation signal that is the target of current in- 
direct detection experiments. A better understanding of such struc- 
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ture within the current concordance ACDM cosmology is thus crit- 
ical for analysing and interpreting all current DM searches. 

We propose a new route to tackle these issues. Rather than 
improving simulations simply by increasing the number of parti- 
cles, we attach additional information to each particle, namely a 
phase-space distortion tensor which allows us to follow the evolu- 
tion of the fine-grained phase- space distribution in the immediate 
neighbourhood of the particle. We introduce the Geodesic Devia- 
tion Equation (GDE) as a general tool for calculating the evolution 
of this distortion along any particle trajectory. The projection from 
phase-space to configuration-space yields the density of the partic- 
ular CDM stream that particle is embedded in and can also identify 
when the particle passes through a caustic. 

This technique makes the fine-grained phase-space structure 
accessible. It enables studies of the phase-space structure of gen- 
eral non-integrable static potentials which reproduce all the results 
previously obtained using frequency analysis methods, identifying 
chaotic regions and finding substructure in regular regions in the 
form of resonances. In addition, it can be used to quantify mix- 
ing rates and to locate caustics. We demonstrate these capabilities 
for the complex phase-space structure of the ellipsoidal logarith- 
mic potential with a core. All relevant phase-space regions could 
be identified by solving the GDE along the orbit. We have written 
a code, DaMaFlow, that allows us to carry out such stream density 
analyses for a wide variety of potentials in a very efficient way. 

Stream density evolution is very sensitive to the shape of the 
underlying potential. We demonstrate this by comparing results for 
a realistic CDM halo with radially varying shape to those for a 
spherical halo with similar radial density profile. After 100 orbits 
the predicted stream densities in the inner regions differ by a fac- 
tor of 100. In general we expect the stream densities to decrease 
as (t/torbitai)~^ for regular orbits and even faster for chaotic or- 
bits, rather than as (t/torbitai)~^. the result found for orbits in a 
spherical potential. Scaling to the Milky Way leads us to estimate 
that there should be at least 10^ streams passing through the solar 
system. 

The potentially revolutionary advantage of our approach, and 
our main reason for pursuing it, is that it applies equally well to 
non-symmetric, non-static situations of the kind that generically 
arise in CDM cosmologies. Indeed, it can be implemented in a rel- 
atively straightforward way in current state-of-the-art cosmologi- 
cal N-body codes. We have carried out such an implementation in 
the GADGET code and have presented some tests based on equi- 
librium Hemquist models. The N-body implementation is able to 
conserve 6-D phase-space density to high accuracy along individ- 
ual particle orbits. In addition, it qualitatively reproduces the results 
found in the corresponding smooth potential for the evolution of 
stream density along individual orbits, and it reproduces the statis- 
tical results found for ensembles of orbits to impressive accuracy. 
The identification of caustic passages is particularly robust, show- 
ing very little dependence on N-body parameters such as particle 
number and softening. Thus discreteness effects appear to be well 
imder control, at least for the large N systems studied here. The 
remarkably robust identification of caustic properties makes us op- 
timistic that we will be able to calculate annihilation boost factors 
due to caustics in fully realistic ACDM dark halos. 

In future applications we will use these techniques to address 
mixing and DM detection issues within fully general simulations 
of the ACDM structure formation model. 
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